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ABSTRACT 

We examine changes in the molecular abundances resulting from increased heating due to a self- 
luminous planetary companion embedded within a narrow circumstellar disk gap. Using 3D models 
that include stellar and planetary irradiation, we find that luminous young planets locally heat up the 
parent circumstellar disk by many tens of Kelvin, resulting in efficient thermal desorption of molecular 
species that are otherwise locally frozen out. Furthermore, the heating is deposited over large regions 
of the disk, ±5 AU radially and spanning < 60° azimuthally. From the 3D chemical models, we 
compute rotational line emission models and full ALMA simulations, and find that the chemical 
signatures of the young planet are detectable as chemical asymmetries in ^ lOh observations. HCN 
and its isotopologues are particularly clear tracers of planetary heating for the models considered here, 
and emission from multiple transitions of the same species is detectable, which encodes temperature 
information in addition to possible velocity information from the spectra itself. We find submillimeter 
molecular emission will be a useful tool to study gas giant planet formation in situ, especially beyond 
i? > 10 AU. 

Keywords: accretion, accretion disks — astrochemistry — planets and satellites: detection — 
protoplanet ary disks — stars: pre-main sequence 


1. INTRODUCTION 

Planetary systems form from the accretion disks encir¬ 
cling young stars. The composition of ice, gas and dust 
within the disk sets the initial chemical conditions of the 
planets and also regulates the physical conditions under 
which planets form. The properties of the pre-planetary 
materials can then be compar ed to the coinposit ion of the 
present day solar system (e.g. Oberg et al. 20111 and even 
that of extrasolar plane tary atmospheres (Maahusudhan 
Teske et aljD013). Nonetheless, there is a missing 


2012 


link between these two stages, separated in time by bil¬ 
lions of years. It is essential to observationally capture a 
forming young planet in situ to put together a complete 
chemical (and physical) history of planet formation. 

The upcoming capabilities of the Atacama Large Mil¬ 
limeter Array (ALMA) will provide extremely high sensi¬ 
tivity and spatially resolved observations of disks (up to 
0.007" at 650 GHz, or 0.7 AU at distances oi d = 100 pc). 
At these scales, detailed disk structure will be readily re¬ 
vealed, and observations of the local environment near 
forming protoplanets should be accessible. Young proto- 
Jupiters are expected to be intrinsically hot as they ac¬ 
crete matter through their circumplanetary disk and lib¬ 
erate gravitational potential energy, thereby generating 
substantial accretion luminosity, Lacc- Theoretical mod¬ 
els of early-stage circumplanetary disks find typical “qui¬ 
escent” accretion levels between M = 10~^° Mrr^ vear~^ 
and 10~^ M f7^ year“^ (Ayliffe & Bate 2009 Lubow & 


Martin 20121. Periodically, the circumplanetary disk is 
theorized to undergo ac cretion outbursts similar to those 
seen in FU Ori objects (Hartmann & Kenyon||1985 Zhu 
et al. 20091, where the planet’s accretion rate jumps to 
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M = (1 — 10) X 10 ® M© year ^. Accretion rates can be 
translated into accretion luminosities, 


L^r.r — 


GMpMp 


( 1 ) 


(Pringle||l981 1, where Mp and Mp are the planet’s mass 
and accretioii rate, and i?in is the inner radius of the 
circumplanetary disk, which is something like the ra¬ 
dius of the planet. Assuming a Jupiter-mass planet and 
i?in = 3 Rjup, the quiescent accretion lev els translate 

Lq (see also Montesinos 
the 


et al. 

^ __ 

2015) 

plane 

can 


During an accretion burst, however 


nosity, where the accretion luminosity due to the planet 
can be as high as Lacc ~ 15 L f^ but for less than 0.1% 
of the planet’s fo rmation time ( Lubow fc Martin||2012 ). 


Kraus & Ireland (2012) reported the detection ot a pos- 


sible embedded protoplanet companion in LkCa 15 for 
which they estimated a substantial accretion luminosity 
of Lacc= 10“^ Lq. Planets with similar accretion lu¬ 
minosities will heat the nearby circumstellar dust disk 
and may give rise to detectable signatures with sensitive 
submillimeter continuum observations, though dete cting 


the circumplaneta ry disk itself will be challenging (|Wolf 
fc D’Angelo|[200^ . 


In addition to the thermal effects on the dust, the plan¬ 
etary accretion heating will have a substantial impact 
on the chemical structure in the vicinity of the young 
planet. In the present work, we explore local heating 
due to a single massive protoplanet (a gas giant precur¬ 
sor) on the three-dimensional chemistry of the surround¬ 
ing circumstellar disk. We do not calculate the chem¬ 
istry of the hot, young circumplanetary disk, but focus 
instead on the larger scale effects on the cold, molecular 
disk within which the planet is entrenched. We vary the 
planet’s orbital location from the star for a fixed accre¬ 
tion rate of M = 10“® Mq year“^, where the physical 
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and thermal structure are described in ^ On this model, 
we calculate the time-dependent chemistry in 2D slices 
over the 3D model, focusing the calculations on th e re¬ 
gion near the planet and opposite the planet (^3.1). For 
species that are particularly sensitive to the presence of 
the planet and have strong submillimeter transitions, we 
identify submillimeter emission line tracers and compute 
the emission based on th e thr ee-dimensional structure to 
determine detectability (13.2). In ^we relax some of the 
assumptions of our model and discuss potential caveats 
to our simple approach, and ^summarizes our findings. 

2. PHYSICAL MODEL 
2.1. Axisymmetric Structure 
2.1.1. Density Model 

The background disk density is described by a simple 
power law in surface density, Eg oc R~^ and vertically 
gaussian with a scale height of h cx 12.5 (r/100 AU)^ AU, 
(3 = 1.1. The gas and dust are uniformly mixed with 
a gas-to-dust mass ratio of 100, and the disk gas mass 
contains Mg = 0.01 M©. The dust i s a simple MRN dis- 
tribution in grain size (ogr cx r~^'^' 


Mathis et al. 


19771 


xgr . g,. __ _^ 

with minimum size of 0.005 /rm and maximum aze of 
1 fj,m. The dust com position is assumed to be Draine 
and Lee astrosilicates (Draine & Lee 1984). We fix the 
central star to have a rnass of M* = 0.5 Mq , an effective 
temperature of T^s = 4000 K and i?* = 2.0 R©, char¬ 
acteristic of low mass T Tauri stars. We consider four 
planet locations, dp = 5 AU, 10 AU, 20 AU, and 30 AU. 
For each location, we simply define a vertical gap in the 
disk corresponding to the Hill radius of a Mp = 1 Mjup 
planet as consistent with theoretical models, which find 
massive planets open radi al gaps comparable to the size 


of their Hill sphere (e.g., Lin fc Papaloizou 1986 


Bryden et al.||199^ |Lubow et al.||1999[ |Jang-(Joiiclelf fe| 


1993 


Turner|2012|). The gap radii (half the gap-width) 
planet models at dp = 5 AU, 10 AU, 20 AU, and 


for the 

30 AU 

are 0.4 AU, 0.9 AU, 1.7 AU, and 2.6 AU, respectively. 
The dust density within the gap has a floor-value at the 
gap midplane of pdust = 10“^^ g cm“^ and is optically 
thin to th e pla net’s infrared radiation, which is described 
in Section [2^ The density within the gap is a factor of 
six orders of magnitude lower than the disk; however, we 
explore higher degrees of gap filling by dust in Section]^ 
We emphasize that these models are designed to under¬ 
stand and isolate the local effects of planetary heating, 
and that we do not include variations in the 3D density 
structure near the pla net, either excess fla ring near the 
edges of the gap (e.g., Jang-Condell 2009), or accretion 
streams, both of which will after the vertical and radial 
disk structure and should be examined in future work 
(see also discussion in Section]^. The density structure 
of the bulk disk is shown in Figure [l^. 

2.1.2. High-energy Stellar Irradiation Field 

We calculate the UV and X-ray radiation field from 
the central star for each of the four planet locations, dp, 
from a single 2D slice and assume the high energy stellar 
radiation held is otherwise azimuthally symmetric. The 
procedure for calculatin g the 2D X-ray and FUV hel d are 
the same as described in|Cleeves et al.|(|2013l|2014a|), us¬ 
ing the Monte Carlo radiative transfer code af Betheff fc 
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Figure 1. Circumstellar disk properties shown for the planet at 
dp = 20 AU. Disk properties away from the gap are similar for the 
dp = 5 AU and dp = 10 AU. The central star is located at R=0 
and Z=0 and the plot is made for an azimuthal angle of 0 = 270° 
(opposite the location of planet). 


Bergin (2011b). The FUV radiation is used to estimate 
the gas temperatures (S. Bruderer, in private communi¬ 
cation), which are coupled to the dust temperatures in 
the region where the planet’s impact is the greatest, and 
as such, our results do not depend on the gas temperature 
in the upper layers where the gas and dust temperatures 
are decoupled. The disk is quite opaque to the FUV 
stellar photons (wavelengths between 930 — 2000 A, in¬ 
cluding that of Lyman-a) due to our use of an unsettled 
disk model. A substantial amount of FUV continuum 
and line (Lyman-a) radiation is down-scattered into the 
gap, aided by the UV-exposed, optically thick outer gap 
edge, which tends to deflect photons both into the gap 
and away from the disk. However, the scattered contin¬ 
uum and line photons within the gap do not penetrate 
far into the gap edges and do not have a substantial im¬ 
pact on the chemistry beyond the illuminated wall. The 
integrated FUV field between A = 930 — 2000 A is shown 
in Figure]^. 

The X-ray radiation field is computed using the same 
transfer code as the UV with the c ombined gas and 
dust X- ray absorption cross sections of |Bethell fc Berglii] 
(2011a) and Thompson scattering. Because the X-ray 
photons are stopped at larger column densities much 
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further into the disk, the presence of the gap does not 
have as significant of an effect as it does for the FUV 
(Figure [^). The gap exposes more of the disk to X-ray 
photonsTout the overall X-ray field is much more sensi¬ 
tive to distribution of disk mass rather than the specific 
gap location; these Hill radius-sized gaps are relatively 
narrow compared to the path length of an X-ray pho¬ 
ton at ^ 5 keV. For co s mic ra y ionization, we use the 
model of Cleeves et al. (20131 for the cosmic ray ion¬ 
ization rate at solar maximum for the present day Sun 
(CcR ~ 2 X 10“^® s“^ per H 2 ) as an upper limit to the 
contribution from cosmic rays, consistent with the upper 
limit determined fro m observations of the TW Hya disk 
( Cleeves et al.|[20T5 1. 


2.2. 3D Thermal Structure 

We calculate the thermal structure of the disk assum¬ 
ing two passive radiation sources: the self-luminous cen¬ 
tral star the protoplanet shining due to its accretion lu¬ 
minosity. We assume a circumplanetary accretion rate 
oi M = 10“® Mq year“^ on an Mp = 1 Mjup and 
Rp = 3 Rjup planet, corresponding to an accretion lumi¬ 
nosity of Lace = 5 X 10“"^ Lq. In this simple picture, the 
planet is treated as a spherical, 1500 K blackbody that 
heats its self-carved gap from the inside. The d ust tem¬ 
peratures are calculated using the code TORUS (Harries 
' 20001 [Harries et aI1|2004[ |Kurosa.wa et al.j 200^ ] [Rin^ 
et al.|200 9D assuming radiative equilibrium with the Lucy 
method (|Lucy||1999|). We note that we do not include 
the differential rotation of the circumstellar disk in the 
dust temperature calculation because radiative equilib¬ 
rium at 50 K is quickly attained within minutes to 
hours ( Woitke||l999 1. For a planet at 10 AU, this times- 
pan is negligible compared to the ~ 30 year period of 
the planet’s orbit, thus planetary heating of the disk is 
expected to remain a local phenomena and will not be¬ 
come highly sheared out by differential rotation. The 
thermal structure for the bulk circumstellar disk (away 
from the planet) is shown in Figure for the planet 
at dp = 20 AU. The large-scale physical properties are 
largely similar for all four planet locations in our mod¬ 
els (i.e., the presence of the gap does not significantly 
change the thermal/irradiation structure away from the 
gap). All four models have an outer disk radius of 50 AU 
and an inner disk radius of 0.2 AU. 

The additional heating by the planetary companion 
primarily results in an azimuthally extended (A0 ^ 
±30°) but somewhat thin (i? < 5 AU) swath of ma¬ 
terial along the gap edges, centered on the planet itself. 
The temperature structures in the midplane and for a 
vertical slice centered on the planet’s location are shown 
in Figure The contours highlight the change in tem¬ 
perature due to the presence of the planet. The planet 
increases the temperature by greater than 10% within 
about AR ~ 1.1 AU, 1.6 AU, 2.7 AU, and 4 AU from 
the edge of the gap for the 5, 10, 20, and 30 AU orbital 
radius, respectively. The absolute change in temperature 
at the gap edges closest to the planet (< 0.5 AU from 
the wall) due to the additional heating corresponds to an 
increase from 42 K to 60 K {dp = 5 AU), 35 K to 47 K 
{dp = 10 AU), 27 K to 34 K {dp = 20 AU), and 24 K to 
31 K {dp = 30 AU) for both gas and dust temperatures. 
These substantial ~ 10 — 20 K changes in cold, dense 


molecular gas will substantially alter the local chemistry 
close to the planet. 


2.3. Chemical Model 

From the physical structure outlined above, we can es¬ 
timate molecular abundances as a function of 3D position 
throughout the disk. The disk chemical code used for the 
calculations is presented in Fogel e t al. (20111 and further 
expanded in Cleeves et al. ( 2U14a|). i'he chernical code it¬ 
self is inherently 1±1U, and so to address this limitation 
we extract 2D azimuthal cuts from full 3D model to com¬ 
pute the abundances, and reconstruct the full disk profile 
from the individual 2D calculations. We consider ten 2D 
slices at 3° intervals in azimuth ranging from the planet’s 
location to 30° away from the planet, and assume that 
the planet is symmetric upon reflection. For the rest of 
the disk beyond 30°, we calculate the chemistry based 
upon one 2D slice on the opposite side of the disk from 
the planet, i.e. the “anti-planet” side. Differential ro¬ 
tation is not included in the present models but will be 
explored i n future work (see further discussion below). 
The 1±1D Fogel et al. (20111 code calculates the time- 
dependent abundances relative to hydrogen as a function 
of radius and height from the midplane based on a fixed 
grid of temperature, density, and radiation field condi¬ 
tions. The chemical reaction network is based on the 


OSU gas-phase chemical network ([Smith et al.|2004|), ex- 

tne 


p anded to include grain-s urface chemistry m the method 
of Hasegawa et al. (19921, where the grain-surface reac- 
tants “sweep-out” the surface at a rate related to the 
binding energy and mass of the reactant. We treat the 
hydroge n binding energy in the method of jCleeves et al.j 
(2014b), where we assume that the binding energy for 
desorption processes is the chemisorbed value such that 
H 2 can form even at high temperatures, while the bind¬ 
ing energy adopted for the rate of hydrogenation reac¬ 
tions on the surfaces of cold dust grains is assumed to 
be the lower, physisorbed value, or L(,(H) = 450 K, such 
that the H-atoms are highly mobile across the ice man¬ 
tle. The reaction network has a total of 6292 reactions, 
including both chemical reactions and physical processes 
(i.e., ionization, desorption, etc.) and 697 species. All 
chemical calculations are examined after 1 Myr of chem¬ 
ical evolution. 


3. RESULTS 

In the following section, we present chemical abun¬ 
dance signatures due to the additional heating from a 
planetary companion as calculated over the disk. We 
then simulate the emergent line emission for the fully 3D 
physical and chemical model. Based upon these emis¬ 
sion models we generate ALMA simulations of planets 
embedded in disks. 

3.1. Chemical Abundance Results 

The most important chemical effect from the planet is 
the thermal desorption of molecular species that are oth¬ 
erwise frozen out as ices in the midplane at the orbital 
distance of the planet, which we term “primary trac¬ 
ers.” As discussed in Section [2^ the dust temperature 
rapidly equilibrates to the local irradiation conditions. 
To assess whether the simplification of neglecting differ¬ 
ential rotation is important for the chemical structure. 












































4 


Cleeves, Bergin, and Harries 


135“ 120° 105° 90° 75° 60° 45° 


135° 



O V 6' > 

■o ■& -o ■•S' -o ■•S' 
R (AU) 



O ^ S' 
R(AU) 


135° 120° 105° 90° 75° 60° 45° 



O •S ■<c>'^S''^0'^''^0 


R (AU) 



R (AU) 



R(AU) 



R (AU) 



R(AU) 



120 g 

96 £ 
72 I 


49 & 


80 g 

65 £ 

50 I 

<D 

35 t 
0) 

20 ^ 



49 £ 
38 I 

<V 

26 t 


<D 

15 ^ 


50 g 
41 £ 

32 I 


24 I 

0 

15 ^ 


Figure 2. Thermal structure in the region near the planet. Left column shows the midplane temperature over radius and azimuth (planet 
is located at 90° in all cases). Right column shows the vertical temperature structure centered on the planet. The thin (thick) contour line 
highlights the region of the disk where the planet increases the local temperature by > 10% (> 30%). Top, middle and bottom rows are 
models for planets at 5 AU, 10 AU, 20 AU and 30 AU, respectively. 


we must compare the chemical timescales for adsorp¬ 
tion/sublimation to the orbital time. If the chemical 
timescales are longer than the orbital time, the chemi¬ 
cal effects of the planetary heating within the gap will 
become sheared out over azimuth. In contrast, short 
chemical timescales relative to the planet’s orbital time 
indicate the relevant chemistry is able to adjust quickly, 
and thus “follow” the planet. 

In our model, the midplane at d = 10 AU has a den¬ 
sity of approximately ~ 10^*^ cm“^ and a temper¬ 
ature of ^ 50 — 60 K near the planet. The timescale 
for freeze-out of a molecule is related to the surface area 
of grains per unit volume, or n((Tgr) cm^ cm“^, as well 
as the thermal speed of the molecule in the gas phase, 
ux, such that molecule X freezes out in a character¬ 
istic time tfo(X) = (n(crgr)ux)~^- At gas densities of 
riH = 10^'^ cm“^, a typical grain surface area density 


is n(agr) 


10 


-11 


cm^ cm ° (i.e., 0.1 /rm-sized grains 


at an abundance relative to H-atoms of 6 x 10“^^), an 
H 2 CO molecule (for example) with mass of 30 amu, 
will collide with a grain on average every ^ 0.3 years, 
which is roughly the time for circumstellar disk chem¬ 
istry to reset when not directly heated by the planet. 
The corresponding timescale for thermal evaporation of 
an H 2 CO molec ule (assuming a desorpt ion temperature 
of ~ 2050 K; Garrod fc Herbst||2006 ) on a 50 K dust 
grain based on the Folyani-W igner relation (see Fogelf 
et al.|[2MT ) is 0.02 years (the timescale for molecules to 


evaporate when exposed to the planet’s heating). Both 
of these timescales are sufficiently rapid compared to the 
orbital time around a 1 Mq star for the planets consid¬ 
ered here (about 11, 32, and 89 years) and therefore we 
expect the chemistry of the primary tracers to follow the 
planet and not experience strong azimuthal shear. It is 
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Figure 3. Vertically integrated column densities of select molecular species (labeled in the lower right corner) after 1 Myr of chemical 
evolution for the planet at dp = 5 AU with an accretion luminosity of Lacc = 5 X 10“^ Lq. The column density at the disk azimuthal angle 
centered on the planet (90° in Figure]^ is shown in magenta, where the points correspond to the specific radii calculated in the chemical 
models. The chemistry of the disk near the gap in the absence of the planet is shown in cyan, i.e., at the ‘anti-planet’ side of the disk. 





Figure 4. 


Vertically integrated column density of select species for the planet at dp = 


10 AU. Figure as described in Figure]^ 
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Figure 5. Vertically integrated column density of select species for the planet at dp = 20 AU. Figure as described in Figure [3| 


important to note that the relevant timescales will in¬ 
crease with height; however, in all cases, the chemical 
effect of the planet is limited to a narrow vertical region 
close to the midplane where zjr < 0.1, over which the 
density only decreases by ^ 30%. 

Molecules that form via gas-phase or grain-surface 
chemistry as a direct result of primary species desorption 
are “second-order” chemical effects - i.e., secondary trac¬ 
ers. The timescales for secondary species formed from 
primaries to return to the low-temperature state will de¬ 
pend on the particular formation pathways for the sec¬ 
ondary species and the availability of He+ ions to break 
up newly formed molecules that would otherwise not be 


present without the planet (e.g.. 

Bergin et al. 112014 

Fu- 

ruya & Aikawa 2014 

). In this case, the distribution of 

sec- 


ondary species formed due to the planet are more likely 
to be affected by shear (see discussion below). 

Given the size of the chemical network, we took an un¬ 
biased approach to search for promising chemical tracers 
of planetary heating, both primary and secondary. We 
calculate the vertical column density of every species in 
the network versus radius, and filter out species with low 
column density near the planet, Nx < 10^° cm“^, since 
these will be the most difficult to detect. We then cal¬ 
culate the fractional difference between the column den¬ 
sities, ANx{R), at the azimuthal location of the planet 
and the anti-planet side, and integrate this quantity over 
radius to get an estimate for the total fractional change 
due to the planet. We then sort by the fractional change 
to identify which species are most affected by the ad- 
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Figure 6. Vertically integrated column density of select species for the planet at dp = 30 AU. Figure as described in Figure [3| 


ditional heating. The results of this “blind search” for 
gas-phase tracers are shown in Figures(dp = 5 AU), 0 
(dp = 10 AU),[|(dp = 20 AU), and§(^^ = 30 AU). 

For certain species, the additional neating changes the 
vertical column density at the edges of the gap by many 
orders of magnitude. For most “signpost” molecules, the 
heating from the planet lifts the local dust temperature 
at the midplane from below the particular species des¬ 
orption temperature to above, causing a large enhance¬ 
ment in the column density at the edges of the gap. Two 
species stand out from this trend for the planet at 5 AU, 
where CS and CH 4 instead exhibit large deficits. The as¬ 
sumed binding energies of CS and CH 4 are Eb = 2000 K 
and Eb = 1330 K, respectively. In the case of CS, the 
5 AU planet causes an initial, fast desorption of CS (and 
higher column density), but the heating simultaneously 
increases the abundance of gas-phase oxygen bearing 
molecules. Over ^ 2000 year timescales, the chemical 
network converts the CS to OCS ice, resulting in the net 
chemical deficit plotted in Figure]^ The orbital time at 
5 AU is clearly much shorter than the chemical timescale 
for this process to occur, and so in this case, the conver¬ 
sion from CS to OCS (and the corresponding CS-deficit) 
is over-predicted in our simplified models. To properly 
model the planet’s effects on the CS chemistry in the 
inner disk requires models that include differential veloc¬ 
ity shear to address these second-order time-dependent 
effects. However, for planets located further out in the 
disk, CS shows simple net desorption behavior without 
further chemical processing for the more distant, 10 and 
20 AU planets. The change in chemical behavior arises 
because the chemistry leading to OCS is not as efficient 
radially further out where the gas-phase oxygen is de¬ 
pleted. As a result, CS becomes a “primary” planet 
tracer for the planets located further out. Thus using 
thermal desorption as a search tool for planets is more 
reliable in the more chemically inactive outer disk. 

For CH 4 , this species’ low binding energy {T^i ~ 27 K) 
allows for it to be in the gas phase regardless of the 
planet’s location. The deficit in CH 4 is purely a second- 
order effect, where the CH 4 at 5 AU primarily forms 
via gas-phase channels originating from CH^. The pres¬ 
ence of additional, sublimated gas-phase species that are 
also able to react with CH^ will form other molecules 
besides CH 4 . These additional pathways stymy the for¬ 
mation efficiency of CH 4 and leads to the CH 4 deficit 
in Fig. of about a factor of three in column density. 
This process likewise occurs over much longer timescales 
than the planets’ orbital times and should be studied in 
further detail including disk differential rotation to fully 
characterize the CH 4 chemistry in the presence of the 
planet. 


There is some variation of the radial extent of the re¬ 
gion affected by the planet as seen in the column density 
plots. For example, as seen in Figures!^ andthe 20 AU 
and 30 AU planets exhibit enhanceoHCbl confined to 
a narrow region close the gap edges (< 0.6 AU), while 
CH 4 is enhanced over a much wider region, < 1.5 AU. 
The main factors that determine the chemically affected 
region are the gap size (i.e., a larger ga p reduces the 
amou nt of heating at the wall’s interface; Cleeves et al.| 


20111 and the binding energy of the species in ques- 
CH 4 is more weakly bound to the grain sur- 
1360 K, comp ared to HCN, Eb ~ 1760 K 


tion. 

face, Eb _ 

(Hasegawa & Herbst||1993|, and so the region near the 
planet is only sufficiently hot to evaporate HCN close to 
the gap walls. 

In addition to HCN and CH 4 , NH 3 is also enhanced in 
the presence of the 30 AU planet (see Fig.|^. This behav¬ 
ior is a feature of the relatively low NH 3 oi nding energy 
used in the present models, Eb = 1100 K ( Hasegawa fc| 
Herbst|1993 ), corresponding to a desorption temperature 


ot 'Jd ~ 28 K. Alternatively, NH 3 deposited on a H 2 O ice 
surface has a much higher binding energy, Eb 3200 K 
(jCollings et al.||2004|), or ~ 90 K. Thus if the higher 
desorption ternperature applies, NH 3 evaporation may 
be a more useful tracer of planets in the inner few AU 
of a cool protoplanetary disk, or further out for warmer 
disks like those around Herbig Ae/Be stars. In summary, 
the ideal chemical tracer for identifying planets is funda¬ 
mentally a balance of the appropriate chemical binding 
energies and the disk thermal structure. 

From the 2D chemical model slices we reconstruct the 
3D abundance profiles. We plot the abundance structure 
for particularly promising species for the four planet lo¬ 
cations in Figures and From these plots, 

it is clear that the memical effects extend over a larger 
azimuthal range than radial range for most cases. The 
abundance enhancement furthermore spans the full wall 
height until the abundance profile merges into that set 
by the surface chemistry, as driven by the central star. 
However, the vertical density structure is such that the 
densities are highest closest to the midplane and most 
tenuous near the surface. Thus the midplane heating by 
the planet can have a large effect on the total, integrated 
column density, similar to the direct midp lane illumina- 
tion of transition disks by the central star (Cleeves et al.j 


20111 . 

UT all the species considered here, HCN is a partic¬ 
ularly robust tracer with a simple midplane chemistry, 
such that the thermal effects from the planet are pri¬ 
marily its desorption. The binding energy of HCN in 
our models is set to 1,760 K ( Hasegawa fc Herbst|1993 ), 
i.e, a characteristic desorption temperature ot ~ 44 K . 
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Figure 7. Chemical abundances for species that are sensitive to the planet’s additional heating. The left column shows the abundances 
in the disk midplane, while the right column shows the vertical abundances centered on the planet. Chemical models are shown for the 
planet at dp = 5 AU. 


For the planets beyond 10 AU, HCN should otherwise be 
frozen-out in the midplane for the majority of the disk. 
Even at 5 AU (midplane temperature of ~ 40 K), HCN is 
just beginning to freeze-out in the absence of the planet. 
Thus HCN should be an excellent tracer of planets at or¬ 
bital distance of dp > 5 AU from the central star. Even 
if the particular HCN binding energy is revised, the loca¬ 
tion HCN as a planet-tracer will scale accordingly based 
upon the midplane temperature profile of the disk and 
the molecules’ binding energy to the grain. Furthermore, 
if the HCN binding energy varies with t he properties of 


the substrate (e.g., 

Codings et al. 

2004 

), the azimuthal 

region away from t 

he planet can 

be used as a baseline 


against which the HCN enhancement can be compared 
independent of the specific binding energy assumed. For 
all four models, the enhancement in HCN is typically an 
order of magnitude in column density due to the plane¬ 
tary heating. In the absence of the planet’s effects, the 
baseline HCN column density on the opposite side of the 
disk is Nucn ^ 10^® cm“^ for the planet at 30 AU and 
A^hcn 10^® cm“^ for the planet at 5 AU. The column 
of HCN in the absence of the planet arises almost entirely 
from the HCN present at the warm disk surface. 


3.2. Line Emission 

Many of these “signpost” molecular species have strong 
rotational transitions that can be used to observationally 
identify and characterize the local physical conditions 
near the planet. Because of its simple interpretation and 
large column density, we focus on HCN emission as a 
tracer in the present paper, but emphasize that other, 
perhaps stronger tracers may exist, and the particular 
tracer will depend on the luminosity of the central star, 
which will set the thermal structure of the midplane at 
the disk radii probed by ALMA. 


We use the flexibl e line modeling code LIME (Brinch 
& Hogerheijde|2010 ) to simulate the emergent line emis- 
sion from the full disk, oriented face-on at a distance 
of d = 100 pc. We present HCN emission results for the 
planets at dp = 10 AU, 20 AU and 30 AU as the HCN sig¬ 
nature for the 5 AU planet was not detectable or distin¬ 
guishable from the background disk in any of the models 
considered. Because LIME samples the model’s native 
grid by selecting random points and accepting/rejecting 
them based upon particular criteria (in our case, normal¬ 
ized density and HCN abundance), the code can poten¬ 
tially not adequately sample relatively small scale local 
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Figure 8. Chemical models for the planet at dp = 10 AU. Panels are the same as for Figure]^ 


features in the disk when considering the full 3D vol¬ 
ume of the disk. Arbitrarily increasing the point num¬ 
ber, however, will add more points into the inner, high 
density regions of the disk that are already well-sampled. 
To address this issue, we cast two grids of points, a large 
grid over the full disk range (including the planet), and 
a refined spherical grid centered on the planet that con¬ 
tains 5% of the total number points (^300,000) used to 
sample the model. This insures that a sufficient number 
of points encompass both the disk and the region near the 
planet. We furthermore tested this technique for a model 
without any planet (but still with the refined local grid) 
to confirm that the refined grid does not introduce emis¬ 
sion signatures. Such sub-grids will be useful in modeling 
all types of substructure with LIME. 

Because HCN is fairly optically thick at the surface, we 
model both the HCN and H^^CN isotopologue assuming 
an isotope ratio of = 60. We emphasize that at 

the high densities where the planet is depositing its heat¬ 
ing that selective ph otodissociation effe cts on N 2 relating 
to HCN chemistry (Heays et al. 2014) are not expected 
to play a major role! Tor the emission calculations, we 
consider the J = 4—3 and J = 8—7 rotational transitions 
for both species, which we identify as likely strong transi¬ 
tions accessible with ALMA using preliminary RADEXj^ 
models (RADEX is a statistical equilibrium solver under 
the assumption of the large v elocity gradient approxima¬ 
tion; ]20Rf| ). Using LIME, we calculate 

the line optical depth tor the four transitions considered. 
For both HCN lines, the disk is optically thick (r > 100) 
at all positions at the line center. However, in the line 
wings at velocities 5v > 0.3 km s the circumstel- 


lar disk is optically thin for both transitions, while the 
emission near the planet is thick. Thus these lines will be 
useful for constraining the temperature of the emitting 
circumplanetary medium. The H^^CN J = 4 — 3 line is 
also thick at line center across the disk (r 15 — 50), but 
becomes thin and the lines, and has an optical depth of 
T ^ 1 near the planet. H^^CN J = 8 — 7 is marginally 
optically thick at line center (r ^ 10), and thin at all 
locations in the wings. 

The LIME models represent “perfect” images of the 
emission lines without any thermal noise or beam- 
convolution (besides the inherent limitations of the pixel 
size, which is 0'.'0025 per pixel or 0.25 AU in our mod¬ 
els). To create more realistic emission models, we take 
the LIME output and, using the simulation capabilities of 
CAS^I^we comp ute simulated ALMA o bservations in the 
same method of Cleeves et al. (2014aI. The alma.outl3 


antenna configuration is used for the dp = 10 AU sim¬ 
ulation of both J = 4 — 3 lines and the alma.outlO an¬ 
tenna configuration was used for the J = 7 — 6. For the 
planet at dp = 20 AU (30 AU), alma.outlG (ahna.outl4) 
was used for J = 4 — 3 and alma.out.il (alma.outOO) for 
J = 7—6. Configurations were chosen to provide optimal 
sensitivity and signal-to-noise to detect the planet. We 
add realistic thermal noise with 0.25 mm precipitable wa¬ 
ter vapor, and then reconstruct the image using CLEAN 
applied with natural weighting. The synthesized beam of 
the simulated observations is typically between 071 and 
(y.'2. The results of the simulated line observations (with 
noise) are shown in Figure!^ for the dp = 10 AU planet, 
Figure 12 for the dp = 20 AU planet, and Figure 13 for 


the dp = 30 AU planet. 


^ http://home.strw.leidenuniv.nl/ moldata/radex.html 


^ http://casa.nrao.edu/ 
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Figure 9. Chemical models for the planet at dp = 20 AU. Panels are the same as for Figure]^ 


In most cases, the emission from the planet is diffi¬ 
cult to disentangle from the emission originating from 
the disk. Ho'wever, the planet’s sublimation effects are 
confined to a ±30° azimuthal extent, -which leaves the 
remaining ~ 300° as a counterpoint against which we 
can compare. To better isolate the signature of the 
planet (without making use of a priori knowledge of the 
planet’s location), we average the emission from the en¬ 
tire disk in annuli to estimate the radial line emission 


profile of the circumstellar disk (Figure 11 12 and 13 
middle row). Then, using the average disk profile, we 
subtract off the smooth background component from the 
simulated observation and examine the residual features. 
This technique will only work if the asymmetric feature 
is small and/or weak relative to the disk’s emission, oth¬ 
erwise a median or averaging azimuthally restricted an¬ 


nuli (to mask out the feature) should be used instead. 
Though this system is face-on, the technique would sim¬ 
ilarly work for systems that are not face-on if the in¬ 
clination is known, i.e., by taking a projected annular 
average. 

TT] and [T^ shows the 
lanetsTo 


The bottom row of Figures 


residual emission, where the planets location is circled in 
all of the panels. Without the additional knowledge of 
the planet’s location, the signature in the J = 4 — 3 is 
not strong enough to identify emission from the noise for 
the planet, primarily because of the high opacity of this 
line, which partially hides the deep emission from the 
midplane. The lower opacity for J = 8 — 7 for both iso- 
topologues, however, makes for a substantially stronger 
signal from the planet, even though the noise is substan¬ 
tially higher at these frequencies. The main limitation 































































10 


Cleeves, Bergin, and Harries 



o cP •^5' 

135° 120° 105° 90° 75° 60° 45° 




14 

-2b)NH3 

12 


10 


< 8 


N 6 


4 

, 

2 

1 

0 

1 1 1 


T 


-7.5 

-9.1 

-10.8 I 

1-12.4" 

1-14.0 


135° 


O <9 

120° 105° 90° 75° 60° 45° 


15 20 25 30 35 40 45 



^ 

R (AU) 



25 30 35 

R (AU) 


Figure 10. Chemical models for the planet at dp = 30 AU. Panels are the same as for Figure]^ 


is the extremely high sensitivity required to measure the 
small fractional change in the emission signature from 
the bright circumstellar disk emission. 

Because in our models the HCN abundance is highest 
along the walls of the gap near the planet, the face-on in¬ 
clination provides the largest column density due to the 
planet. However, if a more inclined viewing geometry 
may allow direct imaging of the gap wall depending on 
the overall scale height of the disk. We examined emis¬ 
sion simulations of a disk oriented at 60° with the planet 
at 90° from the projected rotation axis (in the gap cor¬ 
ner), and were unable to recover any emission from the 
planet. If the planet is on the far-edge of the disk such 
that the outer gap wall is more directly viewable, and 
the emission signature may be observable but only for a 
fraction of the orbit. Thus low inclination (face-on) disks 
are favorable for planet searches with chemistry. 

4. FURTHER CONSIDERATIONS 

To simulate the 3D chemical models presented here, 
we have made a number of simplifying assumptions that 
we discuss further here. The assumed accretion rate 
for our models is on the high end of the range for 
quiescent planetary accretion as predicted by models, 
M = 10“^° —10“® Mq year“^, and is assumed to be con¬ 
stant. We found that it was difficult to detect the chemi¬ 
cal signature of the planet for lower circumplanetary disk 
accretion rates, M < 10“® Mq year“^. However, we em¬ 
phasize that the corresponding accretion luminosity of 
the planets in our models, Tacc = 5 x 10“^ Lq is indeed 
similar to that observed for the tentative LkCa 15 planet. 
Furthermore, if young planets undergo frequent accretion 
outbursts (increasing t he planet’s luminosity b y many or¬ 
ders of magnitude, see|Lubow fc Martin|2012| and if the 


cadence of such outbursts is competitive with the chem¬ 
ical timescales, the region affected by the planet may be 
much larger and more readily detectable. 

One significant simplification in the present models is 
the use of a constant gas-to-dust mass ratio of 100. Disks 
are observed to have a deficit of small dust grains in their 
upper layers, attributed to the settling of grains from the 
surface into the midplane (Furlan et al. 20061. One of 
the overall effects of settling is to create a warmer disk, 
where stellar radiation penetrates deeper. If the disk be¬ 
comes too warm, ice sublimation due to the planet may 
be less effective and one would have to i) observe gas- 
phase tracers with relatively higher binding energies to 
achieve the same contrast between the planet and the 
background disk or ii) look for planets at larger radii 
where the circumstellar disk is cooler. The addition of 
dust mass (in particular dust grain surface area) in the 
midplane will also increase the rate of absorption onto 
grains (i.e., decrease the freeze-out time) and as a conse¬ 
quence, this may help reduce the effects of shear by disk 
rotation. On the other hand, if small grains are quickly 
converted to large grains with low surface area-to-mass 
ratios, then the freeze-out time will increase, and there 
may be a lag in the removal of planet-desorbed molecules 
from the gas phase thus spreading out the thermal effect 
of the planet over a larger azimuthal area (see Section™. 
We have also used a simplified de nsity model in the 


present calculations. Recent work by Dong et al. 
has shown that the dust density distribution wit. 


(20141 
iin the 


gap is different for small (sub-micron) and large (mm) 
sized grains. Furthermore, the small grains readily fill 
the gap while the gaps are clear of the large grain pop¬ 
ulation. In that work, the authors find a dust density 
within the gap of pdust = 10“g cm“^, which is four 
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Figure 11. Select emission lines for the planet at dp = 10 AU. The top row shows lOh simulated observations of the indicated species as 
velocity integrated line emission. The middle row plots the annular averaged emission profile, that is representative of the background disk. 
The bottom row shows the residual between the full disk simulation (top) and the averaged profile (middle), highlighting the asymmetric 
emission. Contours for the residual plots are Scr and 5cr (where applicable). The magenta circle highlights the location of the planet. 


orders of magnitude higher than our simple, relatively 
empty gap model. To explore how the addition of dust 
within the gap affects our results, we have calculated 
models for the planet at 10 AU where we have increased 
the minimum dust density within the gap by four and 
five orders of magnitude to simulate filled-gaps. The en¬ 
hanced dust models have ga p densities that are similar to 
and above that predicted in Dong et al. (2014), and the 
dust density drop within the gap is a factor ot 10 and 100 
below a correspondingly gapless disk. The results of this 
modeling is shown in Figure |14[ The most remarkable 
feature of these models is that the thermal structure of 
the protoplanetary disk near the young planet is similar 
regardless of the amount of dust within the gap. More¬ 
over, the surrounding protoplanetary disk’s temperature 
slightly increases with increasing dust in the gap. The 
reason for this behavior is that some fraction of the ther¬ 
mal radiation emitted from the protoplanet is able to 
more freely escape in the vertical direction for the low 
density gap case and is accordingly lost to space. In¬ 
creasing the dust within the gap naturally traps more 
of the planet’s heating, which is then reradiated isotrop¬ 
ically, sending more heat into the protoplanetary disk 
that would otherwise be lost. We find a small increase 


in the HCN column corresponding to the slight increase 
in protoplanetary disk heating by this effect. The overall 
chemical structure is relatively unchanged. 

Another structural simplification made in this work 
was the assumption of axisymmetry. Massive planets will 
gravitationally interact with the circumstellar disk, gen¬ 
erating density waves that can form spiral arms, chang¬ 
ing the disk’s radial and azimuthal d ensity structure (see 
the review of Kley & Nelson 2012). Planets can also 
excite vortices, creating large azimuthal pressure traps, 
observationally seen as asy mmetries in the dust distri - 


bution (as seen in IRS 48; van der Marel et al. 2013|). 
Furthermore, the presence of the gap exposes tne top 
of the outer gap edge to direct stellar irradiation, which 
enhances heating near th e gap at the disk surface (Jang- 
Condell & Turner 2012). While the latter effect is az- 
imuthally symmetric, tne heating can be locally exac¬ 
erbated by the presence of a massive planet perturb¬ 
ing the vertical disk structure, creating local hot/cold 
spots ne a,r the stellar-irradia ted surface at the planet’s 
location (Jang-Condell 2009). We emphasize that all of 
these planet induced structures/physical effects can hap¬ 
pen in concert with planet-heating and ice-desorption in 
the midplane. The presence of spiral arms can be po- 
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Figure 12. Select emission lines for the planet at dp = 20 AU. Panels are the same as for Figure [TT| 


tentially confirmed as planetary in origin if associated 
with a warm accreting source. Furthermore, the stellar 
gap-edge heating occurs primarily at the surface of the 
disk, away from the planet’s primary effects near the mid¬ 
plane, such that the use of spectral lines with high critical 
densities may help mitigate observational confusion. Fi¬ 
nally, the concentration of dust in pressure bumps can 
alter the chemical timescales such as the freeze-out time, 
though it should be noted that the small grains (which 
dominate the freeze-out given their high surface area to 
mass ratio) are still azimuthally distributed in IRS 48, 
providing a substrate for freeze-out. Thus ideally one 
would find evidence of planets through multiple signa¬ 
tures, both chemical, physical, and thermal (e.g., b y de¬ 
tectin g circumplanetary disk continuum emission; Zhu 

Finally, we have assumed very simple adsorp¬ 
tion/desorption physics, with a single HCN binding en¬ 
ergy. In reality, the HCN ice is expected to be mixed 
with less volatile ices, such as H 2 O or CH 3 OH. Trap¬ 
ping can thus increase the desorption tempe rature to 
that of the strongly bound ice, ^ 5000 K (e.g., Sandford 


et al.|198'^|Sandford fc Al lamandolajl990[|Collings et al. 

2003| ICollings et ai.||2004(), potentially reducing the ef¬ 
fects predicted m this work. How trapping mediates des¬ 
orption will essentially depend on how ices were initially 


deposited on the grain surface substrate - simultaneously 
or sequentially, the thickness of the ice, whether the ice is 
amorphous or crystalline, which impacts the d iffusivity 
of the trapped species (e.g., Collings et al.|2003 l, and the 
specific molecule in consideration, and whether it is “CO- 
like” or “H 20 -like.” Ices like HCN are expected to be 
intermediate species, where some fraction directly subli¬ 
mates and some fraction remains trapped, while, for ex- 
amp le, CS, is expecte d to desorb more readily upon heat¬ 
ing (Viti et al. 2004). While incorporating the detailed 
physics of the desorption process is beyond the scope 
of the paper, they should be taken into consideration 
when comparing to the inherently more complex observa¬ 
tions. We note that if HCN is partially trapped (making 
for a weaker but existent signal), the technique of com¬ 
paring the local (planet) emission signature against the 
disk-averaged emission will nonetheless help compensate 
for more complex chemical effects than considered in the 
present work. 


5. CONCLUSIONS 

We have examined the chemical structure of simple toy 
models of a protoplanetary disk heated by two sources, 
the central star and an embedded, luminous young proto- 
Jupiter accreting at M = 10“® Mq year“^. Chem¬ 
ical species with intermediate freeze-out temperatures, 
around ^ 40 K, will be particularly useful tracers of 
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Figure 13. Select emission lines for the planet at dp = 30 AU. Panels are the same as for Figure [TT| 


the additional planetary heating, as the midplane will 
be colder than 40 K outside of i? > 5 AU. A planet 
in a Hill-radius sized gap will heat the gap edges by 
^ 10 — 20 K, which will sublimate certain species and 
greatly enhance the column density of these species lo¬ 
cally. These changes in column density will translate to 
observable signatures in the circumstellar disk emission, 
which will appear as emission line asymmetries. Such 
asymmetries can be highlighted by subtracting off the 
azimuthally averaged line brightness profile of the cir¬ 
cumstellar disk. We found that HCN isotopologues are 
a particularly robust tracer of the heating effects of the 
planet. Additionally, multiple lines of HCN should be 
detectable, which will provide additional physical infor¬ 
mation (temperature) about the planet’s local environ¬ 
ment. In addition to the excitation information afforded 
by molecular line observations, the use of gas tracers (as 
opposed to continuum) will also provide spectral infor¬ 
mation that can be used to look for velocity substructure 
near the planet. This technique is most useful for plan¬ 
ets between R ^ 10 — 30 AU, where the inner disk is 
innately too warm to provide sufficient temperature con¬ 
trast, and planets in the outer disk are expected to carve 
larger gaps (have larger i?Hiu)) and as a direct conse¬ 
quence the planet’s heating of the disk is more diluted at 
larger radii. 


Finally, the primary (sublimated) species can trigger 
additional chemical processing, creating so-called “sec¬ 
ondary” planet tracers. These secondary products may 
be even brighter than the primary products; however, to 
fully characterize them we must take into account the dif¬ 
ferential rotation of the disk, which is beyond the scope 
of the present work, but will be the subject of a follow-up 
paper. If the timescales for secondary product formation 
are a substantial fraction of the orbital time, the local 
effects may be sheared out over a large azimuthal range, 
i.e., as an “arc.” If the relevant chemical timescales are 
much longer than the orbital time, the emission signature 
of the planet will become a ring (or double ring) at the 
gap edges. However, in the simple case of rapid sublima¬ 
tion/condensation, the chemical signature is predicted to 
“follow” the planet, implying that we will be able to re¬ 
observe and confirm that the emission is associated with 
the planet as it traverses its orbit. 
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panels show the azimuthal and vertical temperature structure, re¬ 
spectively. The bottom panel shows the HCN column density at 
the azimuthal location of the planet, where maximum heating is 
deposited. The fiducial model is shown as the blue contours in all 
panels (pdust = 10“^^ g cm~^). Models where we have increased 
the dust density in the gap by four orders of magnitude (cyan) and 
five orders of magnitude (black) are indicated. 
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